home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_5.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.6 KB  |  73 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.5 (Runge-Kutta-Fehlberg Method RKF45).
  9. % Section    9.5,    Runge-Kutta Methods, Page 461
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements Runge-Kutta-Fehlberg method
  13.  
  14. % for solving the initial value problem
  15.  
  16. %     y' = f(t,y)    with    y(a) = ya
  17.  
  18. %    Define and store the function f(t,y)    in the M-file  f.m 
  19. % function z = f(t,y)
  20. % z = 1 + y^2;
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(t,y)');...
  24.            disp('z = 1 + y^2;');...
  25. diary off;
  26.  
  27. f(0,0); % Remark. f.m and rkf45.m are used for Algorithm 9.5
  28.  
  29. pause      % Press any key to continue.
  30.  
  31. clc;
  32. %    Place the endpoints of [a,b] in  a  and  b.
  33.  
  34. %    Place the initial value y(a) in  ya.
  35.  
  36. % Place the tolerance in  toler.
  37.  
  38. %    Place the tentative number of subintervals in  m.
  39.  
  40. a   = 0;
  41. b   = 1.4;
  42. ya  = 0;
  43. toler = 1e-5;
  44. m   = 14;
  45.  
  46. [T,Y] = rkf45('f',a,b,ya,m,toler);
  47. points = [T;Y];
  48.  
  49. pause    % Press any key to continue. 
  50.  
  51. clc;, clg;
  52. c = 0;
  53. d = 6;
  54. axis([a b c d]);...
  55. plot(T,Y,'g');...
  56. hold on;...
  57. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  58. if m<=30,
  59.   plot(T,Y,'or');
  60. end;...
  61. xlabel('t');...
  62. ylabel('y');...
  63. title('Runge-Kutta-Fehlberg solution to y` = f(t,y)');...
  64. grid;...
  65. hold off;...
  66. shg; pause    % Press any key to continue.
  67.  
  68. Mx1 = 'Runge-Kutta-Fehlberg solution to y` = f(t,y).';
  69. Mx2 = '     t(k)               y(k)';
  70. clc,echo off,diary output,...
  71. disp(''),disp(Mx1),...
  72. disp(''),disp(Mx2),disp(points'),diary off,echo on
  73.